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In this paper we consider gauge invariant hnear perturbations of the metric and matter tensors 
describing the self-similar Lemaltre-Tolman-Bondi (timelike dust) spacetime containing a naked 
singularity. We decompose the angular part of the perturbation in terms of spherical harmonics and 
perform a Mellin transform to reduce the perturbation equations to a set of ordinary differential 
equations with singular points. We fix initial data so the perturbation is finite on the axis and the 
past null cone of the singularity, and follow the perturbation modes up to the Cauchy horizon. There 
we argue that certain scalars formed from the modes of the perturbation remain finite, indicating 
linear stability of the Cauchy horizon. 

I. INTRODUCTION 

The naked singularities predicted in certain solutions to Einstein's field equations pose a threat to the validity 
of Penrose's Cosmic Censorship Hypothesis (CCH); indeed, the CCH forbids the existence of naked singularities in 
generic gravitational collapse. Nonetheless, certain counterexamples do exist in which collapse results in a naked 
singularity. The best known example is perhaps the Reissner-Nordstrom solution (in the case \charge\ < \mass\), 
however other instances would include the Kerr spacetime [291] . spacetimes containing colliding plane waves |10i] and 
spacetimes featuring critical collapse [l^l- These naked singularities suggest the possibility of information from the 
singularity escaping to the external universe, resulting in a loss of well-posedness of the field equations. Fortunately, 
the Reissner-Nordstrom solution also provides a paradigm for the possible saviour of the CCH: perturbations in 
the metric and matter tensor grow without bound when the Cauchy horizon is approached, with the Cauchy horizon 
undergoing a 'blue-sheet' instability and becoming singular itself (see Chandrasekhar and Hartle 0]; see also Dafermos 
P, 5] for the Einstein-Maxwell-scalar field case and Poisson and Israel \2I})- Thus the naked singularity is unstable 
under linear perturbations, and these perturbations are essential to give our spacetime the genericity on which the 
CCH depends; so much so that it could well be that naked singularities are an artifact of the high degree of symmetry 
of the spacetimes in which they are typically observed. 

A class of spacetimes with an additional degree of symmetry in which naked singularities are commonly seen to 
occur is the class of self-similar spherically symmetric (SSSS) spacetimes, for example, certain classes of self-similar 
perfect fluid [1^ and dust [l^ solutions, and the self-similar scalar field Q. In previous work, the authors have tested 
the stability of certain members of this class. In [24[ a scalar field was used to model a perturbation and was allowed 
to impinge on the Cauchy horizon of a SSSS spacetime whose matter tensor was unspecified save for satisfying certain 
energy conditions, and in pointwise bounds are found for a scalar wave impinging on the Cauchy horizon of a 
SSSS spacetime. In (25j we considered gauge invariant metric and matter perturbations of the self-similar null dust or 
Vaidya solution. The present paper represents a continuation of this process, in which we consider the perturbations 
of a more realistic and relevant spacetime, the self-similar timelike dust solution of Lemaitre-Tolman-Bondi (LTB). 
We model metric and matter perturbations to take us away from the high degree of symmetry in the background, 
and, after deriving initial conditions on the axis and past null cone of the origin, allow the perturbations to evolve 
up to the Cauchy horizon. There we see that certain scalars built from the modes of the perturbations remain finite, 
indicating that the Cauchy horizon associated with the self-similar LBT spacetime is linearly stable and does not 
display the 'blue-sheet' instability seen in the Reissner-Nordstrom solution. The 'modes' here refer to the coefficients 
of the state variables following a Mellin transform. We will use a commoving spatial coordinate r and a similarity 
variable y that is a time coordinate in the region between the past null cone of the naked singularity and the Cauchy 
horizon. Then the Mellin transform effects G{y,r) — > g{y; s)r^~^. Then the statement above refers to the behaviour 
of the functions g{y; s): when these satisfy conditions on the axis and along the past null cone that correspond to the 
presence of an initially finite perturbation, they remain finite at the Cauchy horizon. Thus we demonstrate that a 
necessary condition for the linear stability of the Cauchy horizon is satisfied. The corresponding sufficient condition 
would entail demonstrating that the inverse Mellin transform of an initially finite perturbation remains finite. We 
discuss this resummation problem below. 

The perturbation formalism of Gerlach and Sengupta 0, Q which we use in this work is very robust in that it can 
be applied to any spherically symmetric background. Moreover, the formalism has been tailored for the longitudinal 
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or Regge- Wheeler gauge which simpUfies the matter perturbation terms. Thus this formahsm has been used by a 
number of authors in order to describe perturbations of sphericaUy symmetric spacetimes, among them perturbations 
of critical behaviour in the massless scalar field by Frolov 0] and Gundlach and Martfn-Garci'a [13] ; perturbations 
of timelike dust solutions by Harada et al [3, [3; ^tnd perturbations of perfect fluids by Gundlach and Martin-Garci'a 
[ill [l^. These analyses (with the exception of Frolov's) primarily rely on numerical evolution of the perturbation 
equations; there is a gap in the literature with regards to analytic or asymptotic solutions to perturbations of these 
spacetimes. 

In broader terms, perturbations of the metric tensor can be thought of as modelling gravitational waves, an impor- 
tant topic in the current scientific community. This formalism has been used for exactly that purpose by numerous 
authors such as Harada et al [H, [H, [l3| , Sarbach and Tiglio HI] , and similar analyses by Nagar and Rezzolla [13] ■ 
Gravitational waves manifest themselves at the quadrupole and above, that is multipole mode number I > 2. There- 
fore in this work we will consider only those modes I > 2. In addition, we restrict our analysis to the even parity 
perturbations as it is in the even sector where the metric and matter perturbations are fully coupled, thus presenting 
a more substantial and interesting model. In the odd sector, the metric and matter perturbations are coupled but 
only insofar as the matter perturbation acts as a source term, and obeys a decoupled equation that fully determines 
the evolution of the matter perturbation. Furthermore, the master equation governing odd parity perturbations takes 
the form of a wave equation with a source term. This source term is completely and explicitly determined in terms 
of initial data, and does not give rise to any divergence. Then the perturbation may be dealt with without recourse 
to a Mellin decomposition using the methods of [22] and it seems clear that no instability arises. In order to restrict 
the length of the present paper, we defer a complete discussion of the odd parity case to a future publication 

The principal finding of this paper is that the Cauchy horizon formed in the collapse of the self-similar Lemaitre- 
Tolman-Bondi spacetime is stable under linear gauge invariant perturbations in the metric and matter tensors, at 
the level of the Mellin modes as outlined above. In the next section we describe the mathematical background to 
the stability analysis, namely we derive the metric and matter tensor for the self-similar timelike dust solution, we 
outline the perturbation formalism of Gerlach and Sengupta, and we describe two important mathematical tools: the 
Mellin integral transform and the generalized Frobenius theorem. In Section III we test the mode stability of the 
LTB spacetime by finding asymptotic limits for the perturbation modes on the axis and past null cone of the origin, 
and under suitable initial conditions allow the perturbation to evolve to the Cauchy horizon and beyond. We use the 
conventions of Wald [1^ and set G = c = 1 . 



II. PRELIMINARIES 



A. The self-similar LTB spacetime 



The Lemaitre-Tolman-Bondi spacetime has been well studied in the literature and we will not derive this solution 
here (but see for example Harada et al. [l6j). we merely give a summary of the main points: 

The LTB solution describes dust particles which move along timelike geodesies in a spherically symmetric spacetime, 
and thus has a matter tensor of the form 

where u^^u^ = —1. We use comoving coordinates t,r with = 5^ and u^V — 0, and let R = R{t,r) denote the 
areal radius. Solving the field equations gives (letting dot and prime denote differentiation w.r.t. t and r resp.) 

(1) 



87ri?2i?'' 

where m denotes the Misner-Sharpe mass, and in the marginally bound case = ^m{r)^tc{r) — i]^, with tc{r) — 



2 



g \Jr^l2Mi. Thus once we have specified m(r) (or alternatively /9(0, r)) we have completely determined all the unknowns. 

From (IT|) we see the density diverges when i? = 0, that is when t = tc{r). This is the curvature singularity known 
as the shell- focusing singularity, and we can interpret the function tjyr) then as the time of arrival of each shell of 
fluid to the singularity. Note there is an additional singularity known as the shell-crossing singularity when R' = 0. 
We will not consider this singularity as one may extend spacetime non-uniquely through the shell crossing singularity, 
see Nolan [2l|. To rule out the occurrence of the shell-crossing singularity we take R' > for all r > 0, see Nolan and 
Mena [11]. 

Thus the line element for marginally bound timelike dust collapse is 

ds^ = -dt^ + R'^dr^ + R^dn"^. 



3 



We, however, are interested in self-similar collapse, and thus we look for a homothetic Killing vector field which 
solves the equation S/a^b + ^b£,a — '^9ab = 0. If = (a(i, r), f3{t, r)) , this returns four equations, 

d = l, R'^l3 = a', PR" + (3'R' ~ R' + aR' ^0, f3R' + aR-R^0. 

From the first equation we can write a = t + Fi{r), for arbitrary Fi. Since /3i?" + (3'R' — {PR')', we may combine the 
third and fourth equations to give a' = 0, and thus we may change the origin of t to set a — t. The second equation 
therefore gives P = F2{r), and we can make a coordinate transformation to set P = r. The remaining equations are 

t{R'J + r{R')' = 0, tR + rR' -R = 0. 

The first of these equations is ^"^daR' = which is solved if and only if R' is a function of a similarity variable, in 
this case y = t/r. Thus if we set R = rG{y), where G is a function of the similarity variable, we have dR/dr = 
G — y{dG/dy), which is solely a function of y. 

Thus the line element for a self-similar spherically symmetric timelike dust will be 

ds^ = -dt^ + {G~yG 'fdr^ + r^G^dVi^, (2) 

where from now on a prime denotes differentiation w.r.t. y. We may use this metric now to generate the Einstein 
tensor and examine the field equations, still using the co-moving coordinates. The rr component of the field equations 
is G"^ + 2GG" = 0. Integrating yields GG'^ — p^, where p is some constant. The tt component then gives 

1 

^~ 8^r-2G2(G-tjG')' 

which is why we chose GG'^ = > 0. Finally integrating this equation and using -R|t=o = r we can solve for G as 

G(2/) = (l-/.y)2/3, (3) 

where /z = — fp. We note that fiat spacetime is recovered by setting /i = 0. 
There is a shell- focussing singularity therefore at y = Since 

dR/dr = (1 - (1 + _ ^y)-^) , 

we see that prior to the formation of the shell focussing-singularity, y < /i^^ 1 — /i?/ > 0, thus dR/dr > 0. This 
rules out the formation of shell-crossing singularities prior to the formation of shcU-focussing singularities. 

The last issue is to examine the causal structure of the spacetime. Radial null geodesies satisfy 

dt _ ^dR 
dr dr ' 

with the plus and the minus describing ingoing and outgoing null geodesies respectively. Since t — yr this equation 
may be rewritten as 

dr r \ dr J 

If there is some y = constant which is a root of the right hand side of this equation, it represents a null geodesic which 
reaches the singularity in the future/past. Thus the Cauchy horizon, y — yc, is given by the first real positive zero of 

G-yG'-y = 0, (4) 
if one exists, and the past null cone of the origin, y = yp, is given by the root of 

G-yG' + y = 0. (5) 
Since G = (1 — iiy)^^"^, we find there is a Cauchy horizon, and therefore a naked singularity, if /i is in the range 

< M < M*, ^* ^ (^^"^ ~ GOVtj w 0.638014. 

Moreover, when /i is in this range, we have the following: there is one past null cone of the origin yp-, there is an 
additional future similarity horizon dX y = ye > yc', as — > /i*, j/e — >■ y^, and as ^ ^ 0, ?/p ^ — 1, ?/c ^ 1 and t/e oo. 

Thus when < /i < /i* , we have a spacetime with the structure given in Figure [T] The scaling origin at which 
the singularity initially forms is the point {t,r) = (0,0). The apparent horizon forms when g°'''VaR^bR — which 

is equivalent to dG/dy ~ 1. This occurs at y = ^ ^1 — (^)^^, that is, before the formation of the shell-focusing 

singularity y — —. 
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i? = 0,r > 




FIG. 1: Conformal diagram for the self-similar LTB admitting a globally naked singularity. There are three similarity horizons 
at which the similarity coordinate y is null: y — i/p the past null cone, y ~ yc shown dashed, and y — y^ shown as a double 
line. We identify y — yc as the Cauchy horizon, and will call y = ye the second future similarity horizon (SFSH). The apparent 
horizon is shown as a bold curve. 



B. Gauge invariant perturbations 

We will use the formalism of Gerlach and Sengupta [1, [§| . This formalism has been well used in the literature and 
so we will only give an outline here for completeness (but see Gundlach and Marti'n-Garcia [14] or the authors [25l | 
for a more detailed description). 

We perform a 2+2 split of spacetime into a manifold spanned by coordinates — {t, t) denoted {M^, qab), crossed 
with unit two spheres spanned by a;° = (0, </>) coordinates and denoted (5^,7af,). A spherically symmetric spacetime 
will therefore have a metric and matter tensor given by 

g^ydx'^dx'' = gAB{x^)dx^dx^ + R^{x^)^abdx''dx^ , 
t^^dx'^dx" = tAB{x^)dx'^dx^ + \t''^R^{x^)-iabdx°-dx^. 

Capital Latin indices will denote coordinates on Al^, lowercase Latin indices will denote coordinates on 5^, and 
Greek indices the 4-dimensional spacetime (i.e. x^ — {x'^,x<')). R is a function on Al^ and gives the areal radius. 
Covariant derivatives on M^M.^ and are respectively denoted 

giiv,\ = 0, gAB\C = 9ab:c = 0, 

and a comma defines a partial derivative. 

We write a non-spherical metric and matter perturbation as 

ff^f = + hf_,u(t, r, 9, 4>), t^i, = i^y + Atf,^{t, r, 9, 4>), 

where from now on an over-tilde denotes background quantities. The spherical harmonics form a basis for functions, 
and from the spherical harmonics we can construct bases for vectors, 

|y, ; 5, = e^%} (6) 

and tensors. 



[y^ab ; Zab = Y^a:b+^l{l + l)Ylab ] } (7) 
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where we have suppressed the mode numbers I, m, X^ab) = 5(^06 + ^ba) is the symmetric part of a tensor, and eab is 
the anti-symmetric pseudo-tensor with respect to such that eab-.c = 0. Using these, we decompose the perturbation 
in terms of scalar, vector and tensor objects defined on A4^, times scalar, vector and tensor bases defined on S^. 
We write the metric and matter perturbation as 



^ ^fhAsY hAYa \ ^^^^ f AtABY _ ^ AtAYa 



I Symm R^KY-fab + GZ^b) i ' I Symm R^At^Y-fab + At^Za 



where, as mentioned previously, we confine our interest to even parity perturbations; that is, those defined using bases 
Y , Y ai Y^ab and Zab- From these, we construct a set of gauge invariant scalars, vectors and tensors given by 



kAB = h-AB - {PA\B + PB\a) 

k = K-2v'^pA 
Tab = AtAB - tABicP^ - i/^'pc\B - ig^'pc\A 



A 



Ta = AtA-t/-pc-RHt\/^)G 
Ti = At' - [pC /r^){RH\l2),c +l{l + l)(t\/4)G 
r2 = Ae-{R'^tj2)G 



■ (metric) (8a) 
(matter) (8b) 



where pa — hA — ^R^G^a, and va = R,a / R- 

We may then recast the perturbation equations entirely in terms of these gauge invariant quantities (see Appendix 
A). Finally we must consider what to measure on the relevant surfaces to test for stability. As explained in [25|], our 
'master' function will be 

<5P_i = \5^a5^>^\''^, (9) 

where 

-^l I kAB, 0^4 = 

with 

Qo = w''w''Y.,ab. 

where n'^, fh^ = r~' w^(9, (p) and m*'^ are a null tetrad of the background and the * represents complex conjugation. 
We note that SP-i is a fully gauge invariant scalar, being both identification and tetrad gauge invariant (see |25i]'). 



SM>o = ^i^e^kAB, 5^,^^h^h^kAB, (10) 



C. The Mellin transform 



The Mellin transform is an integral transform related to the Laplace transform and is particularly useful for equations 
deriving from self-similar spacetimes. It is defined by 

/>oo 

G(y;s)^M[g(x,r)](r^s):^ g(y,r)r'-'dr (11) 

Jo 

with s G C. For this transform to exist, there will be a restriction on the allowed values of s, typically to lie in a strip 
in the complex plane with cti < Re(s) < (72 • The inverse Mellin transform is given by 

-1 /'C+'ioO 

g(y,r)^M-'[G(y;s)]^— r'G(y;s)ds, (12) 

^TTZ J c—ioo 

where c £ M is such that cri < c < (T2. To recover the original function from the Mellin transform, we integrate over 
the vertical contour in the complex plane of s given by Re(s) — c. We emphasize, as this will be crucial later, that 
we do not integrate over all values of s in the interval <7i < Re(s) < cr2, only over the vertical contour defined by a 
specific value of Re(s) in this interval, which we are free to choose. 

The perturbation equations will reduce to systems of ode's in the individual modes of the Mellin transformed 
variables. Resumming the modes to recover the original function is an extremely complicated task and is beyond the 
scope of this paper (although see [2^ for a discussion) , however we will point out that although the finiteness of each 
mode is not a sufficient condition for the finiteness of the resummed original function, it is a necessary condition. Thus 
we will adopt the following minimum stability requirement: for the inverse Mellin transform to exist on a surface we 
must have each mode of the Mellin transformed quantity finite on that surface. 
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D. Extension to the Frobenius theorem 



The theorem of Frobenius is particularly useful in finding power series solutions to ordinary differential equations 
near regular singular points. Consider the following nth order ODE in f{x); there is a regular singular point at x = 
if the ODE is of the form 

+ + . . . + b,,{x)f{x) = 0, (13) 

with each bi analytic at x = 0. We can Taylor expand each hi about a; = 0, and we denote such as expansion as 



m=0 

The so-called indicial exponents (see below) determine the leading behaviour of the series solutions. It is well known 
that when the indicial exponents repeat the solution must contain a logarithmic term, and when they differ by integers 
the solution may or may not contain a logarithmic term (see, for example, [30]). To clarify the structure of the solution 
when the roots differ by integers we give the following theorem due to Littlefield and Desai p^ . 

Theorem 1 {nth order Frobenius theorem) 

Let f(x) solve an ODE of the form (|13p . Then the indicial equation is 

/„(A) = A(A - 1) . . . (A - n + 1) + 6i,oA(A - 1) . . . (A - n + 2) + • • ■ + 6„_i,oA + b„,o, 

whose roots are the indicial exponents. Collect together the indicial exponents which differ by integers into groups, 
and order the elements of each group as 

{Ai, A2, . . . , Aj, . . .} 

such that \i > Ai+i. Then the solution corresponding to Ai is fi{x) = X]m=o ^™'^™^^^ ' '^^'^ linearly independent 
solution corresponding to Xj is 



f,{x) = K, log^"i X J2 + J] A i^, log(^-^) X J2 



m— 



aA('-i) 



(A - X,)A„ 



whe 



Ki = lim 



(A-A,)(^-^) 



An 



di = Xi- Xj e N, ft = : — , Pi = ft = 1. 

t — 1 



III. PERTURBATIONS OF SELF-SIMILAR LTB SPACETIME 



We will consider only modes I >2. We assume the perturbed matter tensor remains that of dust, 

+ At^y = {p + 5p){u^ + 6u^){uy + 5u^). 

If we write the angular part in terms of the spherical harmonics, 5p = qY and 5Uf^ = QaY + CY,a, then for + 5uf^ 
to be a unit, future pointing, timelike geodesic of the perturbed spacetime, as the conservation equation — 

implies must be the case, we must have 5u^ = j^aY + jY^a for some scalar 7. Additional, we must have 

htt = -27,t. (14) 
Using the Regge- Wheeler gauge (in which hA — G ^ — pa), we may write the gauge invariant matter objects as 

Tab = j , Ta^[^q j, T =T =0. 
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Next we calculate the full set of perturbed field equations as given in (jAl|l . We use the equation k^"^ = to remove 
krr = R'^kft, and we use the tt component of (|Ala[) to define g in terms of the other perturbation variables. 

Thus we have a set of five second order, coupled, linear, partial differential equations in the four unknowns 
{ktt, ktr, k,j}, and the two dependent variables t,r. (Had we removed ku with ku ~ ~'^l,t (1141) . these would be 
third order in 7.) In this set of equations we make the change of coordinates 

(t,r)^(2/ = i,r), (15) 

and then perform a Mellin transform over r, reducing the problem to five second order ordinary differential equations 
in y, the similarity variable, and parameterised by s, the transform parameter. The Mellin transforms of our unknowns 
can be written 

ku=r'A{y;s), hr = B{y- s), k = r'K{y;s), j ^ r'+^ H{y; s), 

thus the four unknowns of our set of ODE's are {A, B, K, H}. 

The future pointing ingoing and outgoing radial null geodesic tangents of the background spacctimc in t, r coordi- 
nates are 

= (i?', 1) , - i-R', 1) 

V2R' V2R' 

respectively, since we restrict i?' > to avoid the shell crossing singularity occurring before the shell-focussing 
singularity. The Ai"^ portion (i.e. neglecting the angular part) of the perturbed Weyl scalars then becomes 

After a change of coordinates, and Mellin transform, we may write each mode of these scalars in terms of A and B. 
We define a new variable D = A + B/{G — yC), and the scalars' modes simplify to 

and SP-i = |(55'o55'4|^/^. We must find solutions to the set of ODE's and use them to evaluate these modes on the 
relevant surfaces. 

We can write this set of second order ODE's as a first order linear system 

Y' = M{y)Y (16) 

where a prime denotes differentiation w.r.t. y, and Y — {A, D, K, H)^ . We note that one of the equations in the system 
is H' ~ —A/2, and thus we have recovered since d/dt — ^d/dy. Due to its length, we give the components of 
the matrix AI in Appendix 2. 

Examining the leading order coefficient matrix near the axis reveals that the axis corresponds to an irregular singular 
point of the system (|16p . with multiple zero eigenvalues, and a number of off-diagonal entries in its Jordan normal 
form; all of which make the system methods used in very unattractive. In any case, we anticipate that the system 
methods would break down when eigenvalues of leading matrices differ by integers, suggesting we would at some stage 
need to decouple an equation in one variable, and use its solution as an inhomogeneous term to integrate the other 
equations. We will sketch the decoupling of an equation in H. 

This system can be written as four first order equations, 

hi{A,D,K,H,A') ^0, h2{A,D,K,H,D') ^0, hsiA, D, K, H, K') ^ 0, h4{A,H')^0. 

We solve the first equation for D = fi{A, K, H, A') and substitute this into the other three equations, giving 

h5{A,K,H,A',K',H',A") = 0, hgiA, K, H, A' , K') = 0, hi{A,H') = 0. 

Combining and to remove K' means we can solve for K = f2{A, H, A' , H' , A"), and we are left with two 
equations, 

h7{A, H, A', H', A'\ H", A'") = 0, h4{A, H') = 0. 
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Finally we remove A for a fourth order ODE in H, 

ho{H,H',H",H"',H"") = 0. 
The other variables can be calculated when the solutions for H are found, as 

A^9,{H'), D^g,{H,H',H",H"'), K = g,{H, H\ H" , H'"), 

and thus we can write the scalars 5\I/o,4 in terms of H and its derivatives. 

Having derived the necessary equations, we pause to outline our general strategy for studying linear stability. 
The perturbation equations comprise a first order system of equations in certain perturbation variables X{y, r) from 
which we can extract physically significant quantities. Of principal importance here is SP-i. We seek to impose initial 
and boundary conditions on the perturbations that correspond to the most general, initially finite perturbation that 
satisfies appropriate conditions at the axis. This perturbation is allowed to evolve up to the Cauchy horizon, and we 
then try to determine if the perturbation has remained finite. The question of what is meant by 'finite' is important. 
A minimal condition is that 5P_i be bounded on the past null cone. This however would allow for an infinite energy 
content over the past null cone - or over a space-like surface an arbitrarily short time to the future of the past null 
cone. This leads naturally to the consideration of the L'^ norm of the perturbation. This is finite if and only if (by 
Planchcrcl's theorem) the norm of the Fourier transform of the perturbation is finite. But the Fourier transform 
is related to the Mellin transform by the complex rotation s ^ iz and r = In p. So we are led naturally to consider 
finiteness of the Mellin transform (a weaker condition than finiteness of the norm) as a minimal condition for 
finitcncss of the perturbation. This is the condition on which we will focus: the perturbation X{y,r) will be referred 
to as finite at time yo if the Mellin transform x{y; s) is finite at time y = yo- 



A. Axis. 



We consider first the axis, y = — oo. We make the transformation y = —1/w to put the axis at w = 0, and then the 
transformation w — z'^ io ensure integer exponents in the series expansions about the axis of the coefficients of the 
difl[erential equation. We find ^; = is a regular singular point of the fourth order ODE in H, which we will write as 

^z-'[/i,+O(^)]^=0, hi^Q. (17) 

i=o 

The indicial exponents near z = are 

{-3,2,-4- Z-3s,-3 + Z- 3s}. (18) 

The ambiguity of the value of ,s complicates matters regarding the position of logarithmic terms in the full solution, 
so we will begin by reminding ourselves of the two conditions the solutions must solve: 

1. The solution must exist; that is we must be able to recover the original function from its Mellin transform. Our 
minimum stability requirement for this to hold is that an acceptable solution is one which does not diverge on 
the relevant surface. 

2. The solution must be such that each mode of ^^'o,4 is finite on the axis. 

Consider the indicial exponent —3. Regardless of the values of the other exponents, the corresponding solution will 
contain at least the series Y^'^=o ■ "^^^ is certainly not convergent, it diverges at ^ = (^o 7^ 0). Thus we 

must not consider this solution. 

In examining requirement 2, we expand the coefficients of H and its derivatives in 5^*0,4 around z = Q, and we find 
the dominant term is 

<5^fo_4 ^ r'-'^z'^H'. 

Consider the indicial exponent Z — 3s — 3. The contribution to the general solution corresponding to this eigenvalue is 

00 

(Logarithmic terms) x (Series) + ^ A^^™"'"'"^^"^, 

m=0 
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where the first portion of this solution depends on the other eigenvalues, and may not even be present. Near z = 0, 
we find (55'o,4 ^ z'^^ due to the second term. We would certainly expect these scalars to be finite on the axis for the 
quadrupole and other modes with Z < 6, thus we must rule out this solution. 

Similarly for the indicial exponent —A — I — 3s, we find S^o,4 ^ near z = 0. Thus we must also rule out this 

solution. 

Finally for indicial exponent 2, we see the solution 

oo 

H=J2 ^mZ"'^' (19) 
m=0 

is convergent (near z = 0) V s, and thus satisfies our minimum stability requirement. Further, the scalars S'^qa will 
be finite on the axis for Re{s) > 1/3. Thus we have found a one-parameter family of solutions near the axis. 



B. Past null cone 

The past null cone, y = yp, is the real, negative root of G — yG' + y = 0, where G{y) = (1 — fiy)"^^^. There is only 
one real root (when < < /U*), and it is parameterized by fi. Thus we may write 

G-yG' + y = {y-yp)F{y), F(t/p) ^ 0. 

yp is a very cumbersome surd, and is quite difficult to work with. Instead, we draw out the nature of the coefficients 

of the _ff-cquation by using G'{yp) = {G{yp) + yp)/yp. We find that setting G' = {G + y)/y makes each coefficient 
vanish, except for the coefficient of the highest derivative. Thus we may write the _ff-equation as 

(G - yG' + y)[mo + 0{y - yp)]H^^^ + K + 0{y - yp)]H''^^ + bo + 0{y - yp)]H^^^ 

+ [90 + 0{y - Vp)]H' + [ro + 0{y - yp)]H = 0, 

where mo, no etc. are the first nonzero terms in the series expansions about the past null cone. 
We may write this in canonical form as 

{y ^ yp)^H(^'> + {y- ypf h{y)H^^^ + (y - ypf b2{y)H" + . . . = 0. 

If the series expansions of the bi about the past null cone are denoted hi = X^^o ^»>J ~ 2/?)'' > then the first few terms 
in the expansions of the hi about y = yp are 

bifi = no/{moFp) 6i_i = . . . 

&2,o = 62,1 =Po/(moFp) 62,2 = •■ • f^20) 

hfi = 63,1 = 63,2 = qo/ (mo-Fp) 63,3 = • • • 

64,0 = 64,1 = 64,2 = 64,3 = ro/{moFp) 64^ = . . . 

where Fp = F(yp). Therefore y = j/p is a regular singular point of this ordinary differential equation, and the indicial 
equation for a fourth order ODE is 

A(A - 1)(A - 2)(A - 3) + 6i,oA(A - 1)(A - 2) + 62,oA(A - 1) + 63,oA + 64,0 = 0. 

Thus the indicial exponents are 

{ 0, 1,2, 3 -61,0 = a }. 

To determine what exactly a is, we note 

Fp = F{yp) = lim ~ + ^ = 1 - ypG"{yp), 
y^vp y-yp 



using rHopital's rule, and thus 
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We note that a = s in the Hmit ^ 0. 

We may find the solutions due to these indicial exponents from the analysis in §11. D. Let us consider first the case 
a ^ Z. We group together the indicial exponents as 

{2,1,0}, {a} 

since cr ^ Z. Then, according to Theorem [Jl the general solution has the form 



m+2 



jn—O 



Ki \ogx Y,Ara{y- yp)"+' + ^ - yp) 



m+l 



m=0 



m=0 



Ki log2 X Am{y - VpT^^ + K2P2 logx Y Bmiy " + Y ^^^y - ypy 



m=0 



m=0 



m=0 



rn—0 



Dyn{y - Vp 



(22) 



where the hi are constants of integration and we have used an overbar to distinguish the K coefficients in the second 
and third solution. 

We see the general solution contains three logarithmic terms, each multiplied by a constant. For the fourth order 
ODE in H we are considering here, these constants are (again from Theorem [T|) 



Ki = hm [X^A2\ = - 



^4,l(&3,l + ^4,1) 

-6i,o)(26i,o-: 



K2 — lim 



^4,1 



(2-&i,o)^ 



where we have set Aq — 1. Crucially, since [(|20p] &3^i — 64^1 
we have a general solution 



0, each of these terms vanish, and thus when ct ^ Z, 



V=Vp 



hi 



00 

E 

m=0 



Am{y - yp) 



m+2 



m=0 



Bm{y - yp 



,m+l 



m=0 



Cm(y - yp) 



m=0 



Dm{y - yp) 



(23) 



with each series linearly independent. Our minimum stability requirement for these solutions will be satisfied for 
Re{a) > 0. 

Now we examine the scalars (5\E'o,4 near the past null cone, and we find we can write 



(24) 



5*0 ~ ciH + C2H' + C3H" + ci{y - yp)H^^\ 

with a similar expression for 5"^!^. The scalars are automatically finite on y — yp for the first three series in 
For the fourth series, we find surprisingly that C3cr(cr — 1) + C4cr(cr — l){a — 2) = for both i5*o and 5'^i\ that is the 
coefficient of the leading term, which goes like [y — j/p)'^^^, vanishes exactly. Thus for finite scalars on the past null 
cone due to the fourth solution, we require only Re{cf) > 1. 

Now let us consider cr g Z. Firstly if cr < 0, the minimum stability requirement is not met and we certainly cannot 
recover 7 from H via the inverse Mellin transform; thus we consider cr > 0. Now we note an important point regarding 
the Frobenius method: if two indicial exponents differ by an integer, the solution corresponding to the lowest index 
may contain a logarithmic term; however if two indicial exponents are equal, the second solution must contain a 
logarithmic term. 

If (J = 0, then there will be a solution which has leading term \n{y — yp), which diverges at the past null cone, and 
thus the minimum stability requirement is not satisfied. If cr = 1, the corresponding leading term is (y — yp) ln(y — ?/p), 
which is finite in the limit y Vp- Thus we only consider cr > 0, when cr G Z. 

When calculating the scalars (5^'o,4, we see from ([M)) that if cr = 1, then H' ^ \n{y — yp), and thus we must discount 
cr = 1. Again, when cr = 2, we find H" ^ \n{y — yp), however for cr > 3 we have (S5'o,4 ^ 0(1). 

Thus for the scalars to be finite on the past null cone, we require s to be such that i?e(cr) > 1, with the exception 
of cr = 2. 
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Re{s) 




a = l 



1 



a = 1 



FIG. 2: The lines a = 1 and a = 1 plotted in the Re{s),ii parameter space for < ^ < fi*. 



C. Cauchy Horizon 



The Cauchy horizon, denoted by y — i/c, is the first real root of G — yG' ~ y — where G = (1 — /^y)^^'^ and, as 
described in §11. A, exists and is unique for < /i < /z*. The situation on the Cauchy horizon is very similar to the 
past null cone: we obtain a fourth order ODE in H with y = yc as a regular singular point, use series expansions 
about y — yc oi the coefficients of the differential equation in the form ([^0]) . and find indicial exponents {0, 1,2,ct} 
where 



When (T ^ Z, all the logarithmic terms in the general solution vanish as at the past null cone. The scalar (5^o can 
be written near y — yc as 



with a similar expansion for (5^4. Again, the coefficient of the leading term due to the solution due to the indicial 
exponent a vanishes, and we find the scalars will be finite on the Cauchy horizon iff Re{a) > 1, when ct ^ Z. 

When (T G Z, we find, for the same reasons as at the past null cone, we must rule out ct < 1; when a = 2 the scalars 
diverge like luijj — yc)', and when a > 3 the scalars are finite on the Cauchy horizon. 

Let us consider first the clearer picture, when neither ct or tr are integers. Both a and a are parameterized by s and 
/i, and thus we can plot the line in the Re{s),fi parameter space where a = 1 and a = I. We give this schematically 
in Figure [2] for < ^ < /i*. 

We interpret this plot so: for every if Re{s) is such that the point {fi, Re{s)) is above the fine cr — 1, the scalars 
will be finite on the past null cone. Similarly, if Re{s) is such that the point {Re{s),fj,) is above the line a ~ 1, the 
scalars will be finite on the Cauchy horizon. As the (7 = 1 line is always below the a = 1 line for < /i < /i* , this means 
that all perturbations which are finite on the past null cone at the level of the modes of the Mellin transform will be 
finite on the Cauchy horizon at the same level, when cr, ti ^ Z. It remains to consider the problem of resummation; 
this is discussed below. 

When cr, (T G Z, the picture is a touch more intricate, due to the fact that a — 2 or a ^ 2 will give a divergence in 
the scalars. Consider Figure [31 and let's choose a particular value for /i, fiQ where < //q < The solid portion of 
the line fi = fiQ represents all the allowable values (from the point of view of initial data) of Re{s) for this fiQ, with 
the exception of where the line intersects a = 2. We see that this line must intersect ct = 2 at some point (/io,s*), 
represented by the black dot in Figure [31 

This point represents a precise value of s for which, if we were to perform the inverse Mellin transform over the 
vertical contour in the complex plane of s given by i?e(s) = s*, the perturbation variables thus returned would 
generate finite scalars (5^*0,4 on the past null cone of the origin, but diverging scalars on the Cauchy horizon. However, 
we maintain this is not enough to conclude the Cauchy horizon is unstable, for the following reasons: 




lim a — s. 



(25) 



d^Q ^ ciH + C2H' + C3H" + Ciiy - yc)H^^\ 
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Re(s) 




2 



2 



s 



2 



1 



1 



1 



FIG. 3; The lines a = 1,2 and a — 1,2 plotted in the Re{s),ii parameter space for < /i < /i« . 



1. From our definition of a (I25p . for ct = 2, a real integer, we require s — s* (z M.. Thus there is only a single, 
isolated point in the s complex plane at which s is such that a — 2, and it lies on the real axis. When 
performing the inverse Mellin transform, we must integrate over the contour Re{s) = s* in the complex plane, 
where < s* < ?2, as in Figure |4l Thus the function 7 is recovered as 



A well known theorem in complex analysis (Cauchy's integral theorem), states that we may continuously deform 
the contour of integration if the region thus swept out does not contain any poles. From our solution for H 
when s — s* (and thus a = 2), we see that the integrand has no poles due to the value of y. That the integrand 
has no poles due to the value of s is a technically very difficult question to address full y, a nd is beyond the scope 
of this paper; however, some analysis in this direction was carried out in Section 6 of [23, and there is evidence 
that no poles would be encountered in the general solution for H. 

Thus when performing the inverse Mellin transform we may avoid the single, isolated point which makes the 
scalars diverge. 

2. The diverging mode corresponds to a single isolated point, that is a set of zero measure, in the s plane. This is 
not generic in any sense; to conclude an unstable Cauchy horizon we would be looking for an extended region 
in the s plane in which the modes diverge. 

3. Note that a = 2 lies between a = 1 and a = 2. Thus the value Re{s) = s* means non-integer exponents in the 
solution for H near the past null cone; that is the solution is non-analytic. From the point of view of critical 
collapse, we would restrict our initial data to only consider analytic perturbations, and thus would avoid the 
diverging mode altogether. However assuming analytic initial data is a very strong restriction that we do not 
feel is warranted in the present case. 

For these three reasons we conclude that the Cauchy horizon formed in the collapse of self-similar timelike dust is 
stable at the level of the modes of the Mellin transform under even parity perturbations with I > 2. 



On the second future similarity horizon, denoted y = ye, we find indicial exponents for the fourth order ODE in H 
as {0, 1, 2, a}, where 




D. SFSH 




lim CT = — 1. 



(26) 
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/ 

Re{s) — s* 




Re{s) 










Re(s) = ?i 


Re{s) = <;2 



FIG. 4: Integrating over a contour in the complex plane of s. 

Again we find the scalars go like {y — UeY^^ ■ We may write 

d — 1 — a{fi)s + l3{^). 

Our initial data confined Re{s) > 0, and it is easily found that for < /i < /x*, both a(/i) and f3{fj,) are always 
negative. Thus the scalars (5^o.4 will diverge on the SFSH for all values of Re{s) allowed by initial data (in contrast 
to the Cauchy horizon). 

We conclude that the second future similarity horizon formed in the collapse of self-similar timelike dust is unstable 
at the level of the modes of the Mellin transform under even parity perturbations with I > 2. 

IV. CONCLUSIONS 

We have examined the linear stability of the Cauchy horizon which may form in the collapse of the self-similar 
Lemaitre-Tolman-Bondi spacetime, due to even-parity perturbations of the metric and matter tensors of multipole 
mode I > 2. We have found that the Cauchy horizon is linearly stable at the level of the modes of the Mellin tranform 
of the perturbation variables. However, interestingly, the second future similarity horizon which follows the Cauchy 
horizon is unstable. 

A crucial question then is whether the same result applies to the full perturbation - that is, to the resummed Mellin 
modes. This is a highly nontrivial question. We note two possible approaches. One would be to try to determine 
the asymptotic behaviour of the solutions of the s— parametrised system of ODE's for large values of |s|, with a 
view to showing that the solutions fall-off at a rate that would guarantee existence of the contour integral giving the 
inverse Mellin transform. Another would be to employ the energy methods of [l^l to directly study the evolution 
of the perturbation without recourse to the Mellin transform. This approach is currently begin used by one of us 
(BCN) to study rigorously the even parity perturbations of self-similar Vaidya spacetime. In both cases, there are 
significant technical obstacles. For the asymptotic analysis, one would need global information about how the different 
independent solutions of the s— parametrized ODE's at different singular points are related to one another. For the 
energy methods, the transition from a scalar wave equation to a first order hyperbolic system gives rise to significant 
additional difficulties, principally in determining an appropriate energy functional. It is hoped that developing the 
appropriate 'technology' for Vaidya spacetime (a simpler case) will yield results applicable to the present case. 

However, we maintain that the results derived here are of physical relevance. We have found that a non-trivial 
necessary condition for linear stability is satisfied. Furthermore, and for example, an initially finite perturbation 
constructed from a finite number of Mellin modes will remain finite when it impinges on the Cauchy horizon. 
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Finally, we note that our results here mirror exactly those relating to the stability of the self-similar Vaidya spacetime 
previously studied by the authors, namely that the naked singularity survives the perturbation but only does so for a 
finite time. This adds further weight to the observation of the authors in [2^1 that perhaps a generic feature of naked 
singularities in self-similar spacetimes is the linear stability of 'fan'-type similarity horizons (the Cauchy horizon) and 
instability of 'splash'-type similarity horizons (the SFSH), to use the terminology of Carr and Gundlach [Ij. 
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APPENDIX A: PERTURBATION EQUATIONS 

We give here the full set of perturbation equations for the gauge invariant quantities defined in §11. B. Note we only 
consider multipole modes I > 2 and thus all equations are valid. 

'2v^{kAB\c - kcA\B - kcB\A + 2~gABkcD^'') - 2~gABV^kj,^^c + ~gAB {^^^^^ +\{Gc°+ G,'')+'R}jka'' 
+2{vAk,B +VBk,A +k,A\B ) + 5As(2i^^l^ + 4z;^i;^ - G^^)fccD - 9ab (^2k,^ 1^ +6v''k,c ^^^ fc^ 

^^^-^ + Gc^ + G,^ + 27^^ kAB = -IQttTab, (Ala) 
kf^J"" + nkc"" - ^^y^^^^) - {k.f +2v''k,c +G:^k) 

+ {kcJ^^'' + 2«^fcci3''' + 2(^'^''' + v^v'')kcD) = -167rT\ (Alb) 
k,A -kAc^^ + kc^^A - VAkc^ = -IGttTa, (Ale) 

kA^ = -IGttT^. (Aid) 

Here TZ is the Gaussian curvature of A4^, the manifold spanned by the time and radial coordinates, and thus equals 
half the Ricci scalar oi A4^; also G^i^ is the Einstein tensor of the background spacetime. 

APPENDIX B: COEFFICIENTS OF THE FIRST ORDER LINEAR MATRIX EQUATION 

We give here the coefficients of the first order linear system Y' = MY of Section III: 

Mil ^G^{l + P- 2s2 + 2sG') - G^y{l + f + 2{s^ + S-1) - ^ + s^) G' 
+ (-1 + 8s) G'^) + G2y2 + f + 2s + (-6 + 31 + 3f + is + 2s^) G' 

+ (-14 + 61 + 6f + s ~ 2s^) G'^ + 2 (-2 4- 5s) G'^) + y^G' {-2 + 1 + f+ 

{-A + 1 + f) G' +{-A + l + P) G'' + (-5 + / + G'^ - G''') 

- Gy^ (^-2 + 1 + 1^ + 2 {-2 + 1 + 1^ + s) G' + 3 {-3 + I + + s) G'^+ 

(-14 + 41 + 4/2 - s) G'^ + 4 (-1 + s) G'"^) 

/G (2G^ (1 + s) - 2G (1 + s) yG' + y^G'^^ (G - yG' - y) {G - yG' + y) 
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Mi2 =G'^{l + l^- 2s^) ~ 2G^ {^2 + 2I + 2P + S- 2s^) yG' + y'^G'^ (-4 + / + 

+ {-5 + 1 + G'^) + Gy^G' (-2 {-2 + l + f + s) + (14 -41- 41^ + s) G"^) 
+ G^?/2 (l + l'^ + 2s+ (-14 + 61 + 61^ + s- 2s2) G'^) 

/G (2G'^ (1 + s) - 2G (1 + s) yG" + y^G'^) (G - yG' - y) (G - yG' + y) (yG' - G) 
Mi3 = (G - yG') (2G3sG' + (-2 + Z + y^ (^1 + G'') + G^y (; + ;2 + 2 (-1 + s + s") 
-5sG'^) + Gy^G' (-2 (-2 + Z + _ s) + 3sG'^)) 
/G (2G2 (1 + s) - 2G (1 + s) yG' + y^G'^) (G - yG' -y){G- yG' + y) 
Mi4 = - 2G'^ (g^ (-1 + s) - 2G2 (-2 + s) yG' + y^G' (2 + G'^) + Gy^ (-3 - s+ 
(_4 + .)G'')) 

/G (2G2 (1 + s) - 2G (1 + s) yG' + y^G'^) (G - yG' -y){G- yG' + y) 

M21 = (G - yG') (G - y (1 + G')) (y^ (-2 + Z + Z^ _ 2G') G' (1 + G') + G^ (Z + 

+2s + 2sG') -Gy(^-2 + l + l^ + 2{-2 + l + l^ + s)G' + (-3 + 2s) G'^)) 
/G (G - yG' + y) (^G"" (1 + s) - 2G (1 + s) yG' + y'^G'^) 

M22 = - 2G^ {l + t^- 2s2) - y4 (2 {-4 + 1 + - G') G'^ (1 + G') + 2G^y {l + f + 2s 
+2 (-2 + 2l + 2l'^ + s- 2s^) G' - (1 + s) G'^) + 2G'^y^G' {4-31- 31'^ - 4s 
+ (13 -61- &f -2s + 25^) G' + 2 (1 + s) G'^) 

+ Gy^G'^ (^6Z + 6^2 + 4 (-4 + s) + (-26 + 8Z + 8/^) G' - (3 + 2s) G'^) 
/2G (2G2 (1 + s) - 2G (1 + s) yG' + y^G'') (G - yG') (G - yG' + y) 
M23=(G-yG')' (-2G2s (1 + s - G') - (-2 + Z + Z^) (1 + g') 
+Gy (-2 + Z + Z2-2sG'-3sG'^)) 

/G (G-yG' + y) (2G'' (l + s)-2G (1 + s) yG' + y^G'^) 

M24= -2G'' {-G + yG') (- (G^ (-l + s))+y2G' (2 + G') 
+Gy (-3-s + (-3 + s) G')) 

/G (G-yG' + y) (2G2 (l + s)-2G (1 + s) yG' + y'^G'^') 
M31 = (-(y2 (-2 + / + /2_2G') G' (1 + G')) -G^ (Z + + 2s + 2 sG') 
+Gy (-2 + 1 + 1^ + 2 {-2 + l + l^ + s) G' + {-3 + 2s) G'^)) 

/G (^2 G^ (1 + s) - 2 G (1 + s) y G' + y^ G'^) 
M32 ={G {l + f + 2s)-{l^ + l-4) yG') /g(2G^ (1 + s) - 2 G (1 + s) yG'+y^G'^^ 
M33 = (-2G2sG'+ (-2 + Z + Z2) y'^G' -Gy (-2 + Z + _ 3s G'^)) 

/G (2G2 (1 + s)-2G (1 + s) yG' + y^G'') 
M34= (2G'' (-(G (3 + s)) + 2yG')) /G (2G2 (l + s) _ 2G (1 + s) yG' + y^ G'') 

M4i=-1/2 M42 = M43 = M44 = 0. 
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